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ABSTRACT 


We present a systematic method for constructing boundary conditions (numerical and physical) 
of the required accuracy, for compact (Pade-like) high-order finite-difference schemes for hyperbolic 
systems. First a proper summation-by-parts formula is found for the approximate derivative. A 
“simultaneous approximation term” (SAT) is then introduced to treat the boundary conditions. 
This procedure leads to time-stable schemes even in the system case. An explicit construction of 
the fourth-order compact case is given. Numerical studies are presented to verify the efficacy of the 
approach. 
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Introduction 


Emphasis on the long-time numerical integration of the fluid mechanics equations has increased 
in recent years. As a result, high-order spatially accurate schemes are favored, because of their lower 
phase error. Such schemes, although they are stable in the classical sense (Lax and CI-K-S stability), 
may exhibit a non-physical growth in time. For a fixed time 7’, these schemes converge as the mesh 
size \ /• — > 0. However, from a practical point of view, in order to achieve leasonable ac<ma<v foi 
large 7', meshes much too fine for the computers available in the foreseeable future art' required. 
Since long-time integrations are encountered in present day computations, it is important to devise 
schemes which are not only classically stable but also time-stable. Specifically, they do not allow a 
growth in time tluit is not called for by tho differential ocjmitions. 

To retain the formal accuracy of a high-order scheme, boundary closures must be accomplished 
with accuracies that are at most one order less than the interior scheme [1], For the scalar explicit 
central-differencing case, Kreiss and Scherer [2] have presented a method for constructing a boundary 
condition of accuracy one order less than the inner scheme such that, a generalized sail) mation-bij-parl s 
property of the differential equation is preserved. Strand [3] has used their approach to construct 
in the scalar case, fourth- and sixth-order central-differencing schemes with boundary closures of 
the appropriate order such that the resulting expression for the derivative' satisfies the summation- 
by-parts property. Recent attempts to utilize these boundary closures to numerically solve a 2 x 2 
hyperbolic system have shown that, in certain cases, an unwarranted growth in lime still results. 

In reference [4], the stability characteristic of various compact fourth- and sixth-order spatial 
operators were assessed using the theory of Gustafsson, Kreiss and Sundstrom (G-K-S) [5] for the 
semidiscrete initial-boundary-value-problem (1BVP). This study showed that many of the higher 
order schemes that are G-K-S stable are not time stable. It was concluded that in practical calcula- 
tions, only those schemes which satisfied both definitions of stability were of any uselulness lot long 
time integrations. Of practical importance was a new sixth-order scheme with fifth-order boundary 
conditions which was shown to be G-K-S and time-stable. Recently, however, it has been found that 
most of the high-order schemes that were time-stable in the scalar case, exhibited time divergence 
when applied to a 2 x 2 system. 

In this paper, we outline a systematic procedure for designing time-stable, as well as G-K-S 
stable schemes of high-order accuracy. The new schemes are guaranteed to be time-stable for any 
hyperbolic system (as long as the system has a bounded energy). The first step in this procedure is 
to construct an approximation to the first derivative (internal plus boundary points) that admits a 
summation-by-parts formula. We rely on the work of Strand [3] for high-order explicit formulations. 
For high-order compact schemes, we derive a new methodology for construction of such schemes. 


Appendix I includes an exposition of the methodology, and a detailed example of the fourth-order 
compact central difference scheme with third-order boundary closures. In section 1, we discuss a scalar 
hyperbolic equation. We show that in general a summation-by-parts formula does not guarantee time 
stability. However, we introduce a new procedure for imposing boundary conditions (simultaneous 
approximation term, (SAT)), that solves a linear combination of the boundary conditions and the 
differential equations near the boundary. This technique is an extension of the techniques used in 
reference [6] to stabilize the pseudo-spectral Chebychev collocation method. It is shown that if the 
approximation of the derivative operator admits a summation-by-parts formula then the SAT method 
is stable in the classical sense and is also time-stable. 

In section 2 we discuss the implementation of the SAT method to systems of hyperbolic equations. 
We show that also in the system case, time stability (as well as Lax stability) is assured by having a 
summation-by-parts property for the numerical derivative operator, provided that the SAT method 
is utilized. 

In section 3 we present numerical results that confirm the efficacy of the SAT procedure even in 
the cases where previous attempts could not attain time stability. It is shown that the theoretical 
predictions for the time stability of the SAT method are realized in practice for both the scalar 
hyperbolic case and the 2x2 hyperbolic system. Finally, an optimization of the parameter t (which 
arises in the SAT procedure) is performed, with regard to efficiency and accuracy. 

1. The Scalar Case 


We consider the scalar hyperbolic equation 

du du 

— = A— 0 < x <1 
at ax ~ ~ 

for which there exists the energy rate 

It Jo u2 ( x >*)<fo = V« 2 (M) -« 2 (o,0) 

For positive A, we have the boundary condition 


(1) 


«(M) = g(t) 

We denote by u a vector of the unknowns (uo(£), U\(t), which corresponds to grid points 

X 0 (= 0),Xj, ...X W (=: 1). 


In this work, we deal primarily with compact schemes for the discretization of the spatial operator 
3 x - For a compact spatial operator, the approximation to the first derivative can be written as 

„ du „ 

p -i- = Qu 

ax 

2 


( 2 ) 


where P and Q are (TV + 1) x ( A'" + 1) matrices. We further assume that: 

Assumption I 

(i) Equation (2) is accurate to order m. Specifically, if we denote by v the vector (v(x 0 , t ), ..., v(x N , 
where v(x,t) e C m and = jAx = ^ , and by v x the values of ((|^)o, •••, (|f )^) T then 

Pv x - Q\ = PT e 

where the truncation error T e satisfies 

|Te| = 0(Ax) m 

(ii) The matrix P has a simple structure (preferably tridiagonal) and is easily invertible. 

(iii) There exists a matrix H , and positive constants p i, fi 2 independent of N such that 

Pi / < HP < ^2^ 

specifically, HP is a symmetric positive definite matrix. 

(iv) There exists a matrix G = HQ such that G + G T has only two elements: p 0 ,o and firyv.w- 
In general we require po,o < 0 < gN,N- 

Assumptions 1 and 2 are common to any useful compact scheme. Assumptions 3 and 4 are specific 
to the summation-by-parts requirement for the spatial operator. 

Equation (1) is now semi-discretized using formula (2) to yield 

^ = AP-'Qu (3) 

at 

Note that assumptions 3 and \ from above admit a summation-by-parts formula in the sense that 

-jr = 9ofl u o + 9 n,nu 2 n (4) 

at 

where 

E(t) = j(u(t),HPu(t)) (5) 

In Appendix I we show how to construct a fourth-order compact scheme that satisfy Assumption 
1 and therefore (4). 

Interestingly, equations (4) and (5) were obtained without imposing the boundary conditions. We 
will use the summation-by-parts property defined in equations (4) and (5) to construct a scheme 

3 


that admits a decreasing energy norm when the boundary condition is imposed. Note that the 
way in which the boundary condition is imposed is important for numerical stability. The most 
common procedure of imposing the boundary conditions (A > 0 ), is to use equation (3) to update' 
the unknowns u lh ..u N , followed by overwriting u N = g(t). This procedure accounts for the fact that 
in a general hyperbolic, system the precise location for each boundary condition is not known until 
alter a characteristic decomposition is performed at all boundaries. This procedure (particularly if 
II is a nontrivial matrix), may not yield the estimate (4) with u N replaced by </(/). In short, the 
imposition of certain boundary treatments may ruin the structure of the summation norm, which 
results in a numerical scheme that is not time-stable. 

A simple counter-example is presented which demonstrates the necessity of careful boundary 
implementation. Consider the scalar equation u t = u x with the boundary condition u N = </(/). The 
semi-discretization in the absence of boundary conditions becomes u ( = A u, where A = P~ x Q. As 

described earlier, once the matrix A is formed, the boundary conditions are imposed. This has the 
('fleet, of p re-multiplying the matrix A by the boundary matrix D. Without loss of generality, we use 
the boundary condition g(t) = 0 in this problem; the resulting boundary operator is the matrix 


Os 


1 0 0 
0 1 0 
0 0 0 


For time stability, the resulting matrix A f = D P~ l Q, rather than the matrix A must exhibit a 
summation- by- p ar t norm. 

For simplicity, we discretize the domain into two even intervals, such that the discrete solution 
vec.toi is ( Uu . Uj . u 2 ) . J he boundaiy condition is imposed at 1 / 2 - A first-order discretization that 
satisfies the summation-by-parts energy norm is 


77/48 

( — 19)/ 1 2 

(— 43)/48 ‘ 


' (-25)/ 16 4 

( — 39)/ 1 6 ‘ 

(— 19)/12 

32/3 

( 1 3)/ 1 2 

; Q 3 — 

-4 0 

4 

( — 43)/48 

(-13)/I2 

53/48 


39/16 -4 

25/16 


Not*' that the matrices P and Q satisfy P 3 = P.J and Q :i = -Q? except for r/ 0(0 and In 

this example, the matrix // is the identity matrix. The characteristic equation for the P :i matrix is 
- lf)2A' ! + 2568A 2 — 5026A + 501 = 0. The symmetry of P 3 and the alternating signs of the respective 
terms in the characteristic polynomial guarantee the positive definiteness of P :i . The discretization 
operator A 3 = P 3 l Q :i can be written as 



. 4 , = 


11/1002 ( — 51 ‘ 2)/ 50 1 1013/1002 

(- 55)/334 (- 1 12)/167 279/334 

2059/1002 (— 2500)/ 501 3061/1002 


ation-by- parts energy norm are satisfied by this discretization, and 


All the requirements of t he summation- by- pa 
a. precise energy norm exists in the absence of boundary conditions 

The combined operator A\ = D :i A* becomes 


A = 


11/1002 ( — 512)/501 1013/1002 

(— 55)/ 334 ( — 1 12)/ 167 279/331 

0 0 0 


for which the characteristic polynomial is - 1 002A" -66 1 A 2 + 176A-0. Tin- roots of the rharae o, ,s 
polynomial are A = -0.86317..., anti A = 0.20349..., respectively. The non, erica! so, ut, on wrll R rtn 
i„ Pone as a rest.lt of the eigenvalue it. the right half of the complex plane (RII-I ) ami w.ll not 

time-stable. 

As demonstrated by the previous counter-example, a spatial operator which satisfies the summation- 
bv-parts energy .torn, may not be time-stable. Many of the high-order scl.en.es that sat.sfy the sum 
.nation property are time-stable for the scalar case. A notable exception ,s the s,xth-order <^p ice 
scheme with fifth-order boundary conditions reported ... the work of St.and [. ]. (. « ' PI" 1 '' • 
details of this scheme.) For this sixth-order scheme, time stability can be guarantee, on V . « as 

row and column of the matrices IIP and HQ are removed before matrix inversion and mult.pl, oat, on 

arc performed. 

The underlying reason for the growth in time is the imposition of the boundar y condition opecatm, 

which has an effect on the structure of the norm matrix P ,n D P Q • • I" " ,c * >’ 

destroys the structure of the norm P. In the scalar case, this problem can be dim, natal certa.n 
circumstances. For instance, if the matrix P is a restricted full norm, then D P~ st, II produces a 

useful norm by eliminating the zero element. A restricted orn, is defined where the d, agonal 

is the only nonzero element in the first (or last) row and column of the matrix P (See fit, an, [■ ] • 
A special case of the restricted full norm is the diagonal case, which ,s of some pract.ca in e,es .. 
Ihifortunately, even for cases where P ,s a restricted full norm, stability cannot b. r r , 

case of a hyperbolic, system. An alternative means of imposing boundary cond.t.ons must he fount 

for these rases. 

A, this point, we introduce the SAT methodology for boundary implementation. We show in 
the following text that the SAT method leads not only to stability hut also to „„m«Ub,W.y fa • « 
scalar wave equations, and this property applies to arbitrary hyperhohe systems. 3 he SA „„ > 

5 



involves the indirect imposition of the physical boundary conditions. This is accomplished by adding 
a term to the denvative operator, which is proportional to the difference between the discrete value 
U N and the boundary term g(t). Thus, we propose the discretization, 


du 

“~fa = '^ ( ? u ~ T ^ 9N,nS(u N - g(t )) 


where 


( 6 ) 


s = W-'(0,0,...,0,1) 7 ' (7) 

Contrary to the common Practice of satisfying the boundary condition directly by imposing u N = j(i), 

Tb T ° ' nV “ ves solvl,1 S a derivative equation everywhere, including the boundary points 

e extra term whrch ,s added accounts for the boundary information to within the accuracy of the 
original discretization. Note that the SAT is added no, only to the boundary equation but to other 
points depending on the structure of the vector S, (which is the last column of the matrix H~‘) 

e extra SAT term does not alter the accuracy of the scheme, since the SAT term vanishes upon 
substitution of the analytic solution. 

We now demonstrate that the SAT method yields a Lax stable and time-stable scheme. For the 
tlTbWn ^ MalySiS ' We tal<e S{t) = ^ Pre - multiply «*">«<»■ ( 6 ) ^ « and use equation (7) 


dii 

Hp ~fa = - t ^9n,n(0,0, ..., 0 , 1) t u n 


We now define the energy E(t) as in equation (5) to get 


( 8 ) 


( 9 ) 


dE(t) 2 

dt ~ g°’° u ° 9 n,nu% — rg N<N u 2 N 

With <7o,o < 0 < g N , N , we can immediately state the following theorem. 

Theorem 1.1: 

The SAT method presented in equation (6) is both stable and time-stable if 

T -‘ ( 10 ) 

In addition to proving the stability of the SAT scheme defined in equation (6), we must show that 
he procedure preserves the order of accuracy m of the spatial operator. This is accomplished by a 
direct convergence proof showing that the SAT term indeed preserves the spatial order of accuracy. 

Denote by v the vector (u(x 0 ,t ), ..., u(x N ,t)) T , i.e. the values of the true solution of (1) at the 
grn Points. Combining the accuracy condition found in Assumption / with equation (6) we have 

dv 


P 


fa ~ ~ T ^d(V,rvS [U^M, t) - g(t)} + PT e 

6 


( 11 ) 


Note that u{x N ,t) - g{t) = u(l,<) - g(t) = 0. Now define 

£j(t) = u(xj,t ) - uj(t) 


where Uj(t) solves (6) , to obtain 


P^ = A Qt - r\g N , N St N + PT e 
at 


( 12 ) 


where T e is the truncation error defined in Assumption 1 . We now use the energy estimate presented 
in (9) to obtain 


d ^' HP ^- < (e, H PT e ) 
dt 


and the inequality 


to obtain 


(e, HPT,) < 




(13) 


By assumption I ,the truncation error is of order m, and we get 


y/idHPe) < 0(AxY 


which proves the convergence of the scheme. 

In conclusion, a precise means is now available for the scalar case to impose boundary condi- 
tions that are guaranteed to be time stable, and that preserve the formal accuracy of the original 
discretization. 


2. The Hyperbolic System 

In this section, we explain how to use the SAT method for systems of hyperbolic equations and 
show that the resulting scheme satisfies an energy estimate similar to the one obtained for the scalar 
differential equation. First the system of differential equations is described. 

Let u 1 and u 11 be the two function-valued vectors 

u 1 = (u (1) (x,t),...,u (fc) (x,f)) 
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( 14 ) 



u 11 = (u (t+1) , .... u( r >(;E, /)) 
that solve the system of differential equations 

(ju 1 , du 1 

— = A — — 
at ox 


15 ) 


dii 11 kU thP 

di dx 

where A 1 and A 11 are diagonal matrices of the form 


A 1 = ding( Ai,...,A fc ) 


(1C) 


A /7 = diag{ X k+u ...,X r ) 

In order to impose the boundary conditions we assume that 


A] > A 2 !> ... > A* > 0 > A* +1 > ... > A r 


For this case, a well-posed set of boundary conditions is given by 


U *(M) - /?u II (i,t) + g I (t) 


u II (0,/) = Lu I (0,f) + g II (t) 

Where 

g 1 ^) = (g (1) (t), g (k >(t)) 

and, 

g n (t) = (g (k+1) (t),..., g (r) (t)) 

In equation (17), the matrix R has k rows and r - k columns, while the matrix L has r - k rows 
and k columns. Without loss ol generality, for the stability analysis we will assume that both g*(t) 
and g n (t) vanish. 

Equation (17) is well-posed for any L and R. However, to guarantee no growth in time some 
conditions must be imposed on the matrices L and R. These conditions are 
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Condition I: 


|i||fl| < 1 (18) 

where the matrix norm is defined by 

\A\ = p(A T A)± 

and p(A) is the spectral radius of A. 

The Continuous case 


It is instructive to establish and prove an energy estimate for the continuous hyperbolic system 
although such a proof is well known. The same basic steps that are used in the continuous proof will 
be used later in the text to prove the energy estimate resulting from the semi-discrete hyperbolic 
system. 

Condition I is a sufficient condition for the solution of equation (15) to be bounded in time. In 
fact one can state 


Theorem 2.1: 


Let u l (x,t) and u **(x,t) be the solution of equation (15) with the boundary conditions (17). 
Recall that we take g / = g 11 = 0. Suppose that L and R in equation (17) satisfy Condition I. Define 
an inner product 

(u; ? t?)= f w(x,t)v(X)t)dx (19) 

Jo 

and an energy function E(t) 


E(t) 


e a . 

i = 1 A * 


w (u ,> 


,<0 


) + 


E — ( 
^ | a/ 


u 


(') „(*)> 


<=/t+i 


( 20 ) 


then the time rate of the energy function satisfies 


dE 

dt 


< 0 


(21 


) 


Proof 


We start by differentiating equation (19) with respect to t. to obtain 


dt 


2 /' 

Jo 


dx 


9 



Using equation (15) we obtain 


d(u^\u^) 


dt 


= 2 / XiU^ dx 

Jo 


so that 


d(u^\ u^) 


dt 


= A,( u ' i '(l,() 2 -u"»(0,() 2 ) 


( 22 ) 


Differentiating equation (20) and substituting equation (22) we obtain the energy rate for the system 
as 


f =£|£|(uM(1,0 2 — <’W)- t, |A|(u'->(1,<) 2 -k«(0,() 2 ) 

ai 2=1 i = A: H- 1 


(23) 


relating the time rate of change of the energy function to the energy that crosses the boundaries. 
Note the change of sign in the second term which results from the negative sign of the eigenvalues 
A i for k < i. We must now quantify the magnitude of the boundary terms in equation (23). 

Replacing the sums in equation (23) with the vector operations 

k 

yv i) (M) 2 = u I (i,t) T u I (i,t) 

2 = 1 


E u (l) (0,t) 2 = u II (0,f) T u II (0,i) 

2 = & + l 

we can now make use of the boundary conditions in equation (17) to obtain 

u X (l,2) u I (l,2) = u II (l,f) i?" r /2u II (l,2) 


(24) 


u II (0, t) T u II (0, t) = u x (0, f) T L T Lu I (0, t) 

Substituting the equations (24) and (25) into (23) we obtain 

= u ll (l,t) T {R T R\L\ - |R|}u n (l,t) + uI (0, t) T {L T L\R\ - |L|}u x (0,t) 


Because Condition 1 ensures that 


and 


R t R\L\ - \R\I < 0 

L t L\R\ - I\L\ < 0 

10 


(25) 


( 26 ) 


equation (21) is established. Therefore the continuous energy function E{t) is bounded in time. This 
completes the proof of Theorem (2.1). 

The Semi-discrete Case 


We are ready now to discuss the implementation of the SAT technique for the system in equation 
(15) with the boundary condition given in equation (17). As in section 1 we denote by u a \ectoi 
unknowns (4°, u? , ...u$) T which correspond to the grid points *„(= 0), x u ...x N (= 1). We assume 
that we have matrices P , Q and H such that the scalar case admits a summation-by-parts energy 
norm given in section 1. The SAT discretization of equations (15) - (17) is chosen as 



XiQu 1 — — (Ru )jv 




1 < i < k 


(27) 



where r is a stabilizing 
of the vectors 


. - soArSW^ 0 - (to 1 )?’ - »“') k+l<i<r 

factor to be determined later. As in the scalar case, we choose S w to be one 


5(0 = //-*( 0,0, ...,0, l) r \<i<k 


(28) 


S {,) = H~ *(1,0,...,0,0) T k+\<i<r 

We recall from the scalar case that HP is symmetric positive definite and HQ is skew symmetric 
except for the terms g 0 ,o = {HQ)o,o < 0 and g N , N = ( HQ)n,n > 0- Thus equation (2f) is well 

defined. 

Before proving the stability (and time stability) of the SAT method in equation (27), we would like 
to comment on the role of the matrix H. Explicit knowledge of H is required for the implementation 
of the SAT method, specifically the knowledge of 30,0 and g NtN as well as the vectors S {,) are needed 
to implement equation (27). Thus H is not only a theoretical tool (as in reference [2]) but is also of 

practical importance. 

We are now ready for the stability proof of the SAT method in equation (27). The proof is 
analogous to Theorem 2.1 with the continuous integrals replaced by discrete sums. The scalar 
product is defined, analogous to equation (19), as 

1=0 

11 


( 29 ) 


A different scalar product to be used later, analogous to equations (24), is 



k 


= £ 


U 77l U 77l 


i= 1 


(30) 

! = t'+l 

for m = 0, iV. 

Theorem (2.2) 

Let the SAT method defined by equation (27) satisfy Assumption 1, for the discretization of the 
hyperbolic system defined in equation (15) with boundary conditions (17), (with g r (t) = g n (t) = 0). 
Then the discretization is both stable and time-stable provided that 


i«uii 

Moreover, let the discrete energy be defined 


2 - 2^1 - <t< 2 + 2^/1 - |«|| i | 




as 


MO = £ V (u\ HPu*) + ± M(u'>Pu i ) 

*= 1 ^ 7 — T-t-1 \K\ 


l=k + 1 l^« 

where the scalar product (u 1 , u ! ) is defined in equation (29). Then 

dEft(t) 


dt 


< 0 


(31) 


(32) 


Proof 


As in theorem 2.1 we differentiate the scalar product (u\ f!Pu') and use equation 


^(u\ HPn') = A,(u \HQu*) - g N , N X t T(u $ - (/?u n )5?)(u\ HS&) 


ion (27) to obtain 


1 < i < k 


(33) 


dt 


V, HPW) = Ai(u‘, HQu‘) - s „,„A,t( u !,‘ i - (fiu'lPKu’, HS®) *+l<i< 


We now use the definition of 5<‘ ' from equation (28) and the properties of HQ from Assumption 
to obtain 

d 


dt 


(u\ HPu') = 


12 


</o,oAi(u^ !) ) 2 + A ,</w,n(u{v ) 2 - A igN,NT(u$) 2 + \ x9n , n tu"(Ru")", \ <i < k 


( 34 ) 


(It 


— < 70 , 0 1 Ai | ( W q ^ ) 2 — \\i\gN,N(u\J) 2 + lAil^o.o^Wo *) 2 — |A,|(/o,o7'Wq ^(Lu 1 )), \ k + 1 < t < /' 

Note that in equation (34) we used the fact that the A; are negative for k + 1 < i < r. We must 
now quantify the magnitude of the boundary terms in equations (34). If the sums in equations (34) 
are replaced with the vector operations defined in equations (30) we get an estimate for the discrete 
energy rate — 


,(*)\2 


(i i\HPu l ) = 
(i)Y2 


.('■),/ „IW) 


dEw{t) ... r I Ti , ir i /i _\r..I ..Ii i \ r \ „ r..I d..Hi 


dt 


= |L|<jfo,o[« , u A ] 0 + \L\g NtN (l - + kkv„v[u ,Ru 


JV 


+ |/?|(r - iWtu 11 ,!. 11 ^ ~ |f?kv.w[u II ,u II ] w - </o,o|^|r[u n , Lu^o 
Substituting the estimates 

[u 1 , Ru ll ] N < |u I | N |/?||u II | Ar 
[u 11 , Lu I ] 0 < ju II io|/.||u I |o 

where 


lu^m = 

into equation (35), and collecting like terms yields 


dEjv(t) 

dt 


< - Olu'lJ, - r|i||fi||u I | w |u“| A , + |ft||u“&} 


I| I.JIl 


. 1 1 1 2 


(35) 


+go,o{mr ~ l)|u u |^ - x|Z,||/?||u 1 | 0 |« 11 |o + kHu 1 ! 2 } 

For to be negative we require each curly bracket to be positive. Thus we need 


Mr - l)|u'|J, - T|i||fi||u 1 |„|u ,, |« + |ft||u‘‘|J, > 0 

and also 

\R\(t - l)|u"g - T|i,||fl||u I |„|u "| 0 + |i||u'g > 0 

Both inequalities are satisfied if 

\R\\L\t 2 < 4 (r - 1) 

and this is equivalent to equation (31). Thus, the proof is established. 
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Ii 1..II1 


.II 12 


(36) 



3. Results 


Conventional Boundary Conditions 

Three high-order spatial discretizations (two explicit and one compact) are the focus of the results 
section: the fourth-order explicit scheme with third-order boundary conditions, the fourth-order 
compact scheme with third-order boundary conditions, and the sixth-order explicit scheme with 
fifth-order boundary conditions. All satisfy the summation-by-parts requirement in the absence of 
physical boundary conditions. The fourth-order explicit scheme is reported elsewhere (see [3] or [7] 
for specific details) and will not be derived here. The fourth-order compact scheme is new, and 
a systematic procedure for deriving both it and other compact high-order schemes is presented in 
Appendix I. The sixth-order explicit scheme was first reported in reference [3], but is also included 
in Appendix II. 


First we demonstrate that all three schemes behave in accordance with their respective order 
properties. We then comment with regard to the sixth-order explicit scheme, that satisfying the 
summation-by-parts energy norm is not sufficient for time stability. 


The model problem used to test the three schemes is the scalar hyperbolic equation 


Ou 3u 
dt dx 


= 0, 0 < x < 1, t>0 


(37) 


u(0,t) — sin27r(—/), / > 0 


(38) 


u(:r,0) = sin27r(x), 0 < x < 1, (39) 

The exact solution is 

u(x, t) — sin 27r(x — £), 0 < x < 1, t>0 (40) 

For all calculations, the time discretization used was a fourth-order Runge-Kutta (R-K) method 
with the time step small enough such that the temporal errors are much smaller than the spatial 
truncation error. In all cases, the boundary condition was implemented at the end of each R-K stage 
by overwriting the value of the solution at the boundary point. 

Table I shows a grid refinement study performed on equation (37) for all three spatial discretiza- 
tions. Both the absolute (log L 2 ) error at a fixed time T and the convergence rate between two 
successive grid densities are plotted. 
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\ 



(fourth 

explicit) 

(fourth 

compact) 

(sixth 

explicit) 

Grid 

log L 2 

Rate 

log L 2 

Rate 

log L 2 

Rate 

21 

-0.501 


-1.418 


1.379 


31 

-2.080 

8.96 

-2.133 

4.06 

1.048 

1.88 

41 

-2.607 

4.22 

-2.627 

3.95 

0.137 

7.29 

61 

-3.329 

4.10 

-3.316 

3.91 

-1.302 

8.17 

81 

-3.832 

4.03 

-3.806 

3.92 

-1.798 

3.96 


Table I: Grid convergence of three high-order schemes on u t + u x = 0. 


This refinement study suggests that all three schemes are Lax stable (the exact solution is approached 
at a fixed time T as mesh is refined) and grid converge consistent with each respective theoretical 
rate. The convergence rates for both of the fourth-order schemes asymptote to the theoretical value 
of 4. The convergence rate of the sixth-order explicit scheme is sporadic but is approximately 6 
(5.28 for the interval between 21 and 81 points). This spurious behavior results from the exponential 
divergence of the solution for long times T. At T = 70, the absolute error of the two fourth-order 
schemes is comparable; however, that of the sixth-order scheme is two to three orders of magnitude 
larger. 

These numerical results indicate that the two fourth-order schemes are time-stable; the sixth- 
order scheme is not. Nothing in the definition of Lax stability precludes exponential divergence of 
the solution for long times T as long as the divergence rate is bounded independently of the grid 
used. (See reference [4].) The numerical divergence of the solution results from a spatial operator 
matrix which has an eigenvalue with a positive real part (an RH-P eigenvalue). For long times 7\ 
the solution is dominated by this eigenvalue. 

To quantify this assertion, a comparison is presented between the numerically observed divergence 
rate, and a theoretical prediction from eigenvalue analysis. By assuming that the numerical error can 
be represented as ejv(t) = e^(0)e ftN< , a growth rate a^ is determined. Similarly, an effective growth 
rate as defined by = \G max (At)\ M , is calculated from an eigenvalue determination. (See 

reference [4] for details). Table II shows a comparison of the observed growth rate of the sixth-order 
explicit scheme with the rate predicted from an eigenvalue determination. 


Grid 

& Numerical 

a (S maj: ) 

21 

0.1672 

0.1673 

31 

0.1879 

0.1886 

41 

0.1880 

0.1879 

61 

0.1659 

0.1746 

81 

0.1785 

0.1808 
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Table II: Numerical vs. Theoretical Growth Rate for the sixth-order explicit. 


The agreement is very good, with a slight discrepancy in the comparison on the 61 and 81 grid-point 
cases. 


The time-divergence seen in the sixth-order scheme is the same as that predicted in the counter- 
example presented in section 1. Specifically, numerical time stability is not guaranteed by a dis- 
cretization which satisfies a summation-by-parts property. Very specific boundary treatments must 
be used to guarantee time stability. 

SAT Boundary Conditions (Scalar) 

The SAT method for treating the boundary conditions guarantees time stability for the hyperbolic 
system. This method relies on a spatial operator that satisfies the summation-by-parts energy norm 
for the scalar case and on very specific boundary treatments to ensure time stability. 

We begin by showing that the procedure does not destroy the formal accuracy of the spatial 
discretization. This result was proven in section 1 for the scalar case. Tables III. a and Ill.b show a 
grid convergence study of the SAT method on the scalar wave equation defined by equations (37), 
(38) and (39). Fourth-order R-K time advancement is used for all runs with a time step such that 
no appreciable temporal error accumulates. All calculations are run to time T = 10. In all cases, 
the calculations remained bounded on all grids (and CFL’s less than CFL ma;r ) for times as large as 
T - 1000, which indicates time stability. This result is consistent with the results from eigenvalue 
determinations in which no RH-P eigenvalues were found. 


r = 1 

(fourth 

explicit) 

(fourth 

compact) 

(sixth 

explicit) 

Grid 

log L -2 

Rate 

log L 2 

Rate 

log L 2 

Rate 

21 

-1.2289 


-1.4005 


-2.5750 


31 

-2.0878 

4.88 

-2.0479 

3.67 

-3.8300 

7.13 

41 

-2.5784 

3.93 

-2.5096 

3.70 

-4.6500 

6.56 

61 

-3.2211 

3.65 

-3.1689 

3.74 

-5.7880 

6.46 

81 

-3.6806 

3.68 

-3.6464 

3.82 

-6.6056 

6.54 


Table III. a. Absolute error (log/^) and convergence exponent with SAT parameter r = 1, for the 
fourth explicit, fourth compact and sixth explicit spatial discretizations. 
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T = 2 

(fourth 

explicit) 

(fourth 

compact) 

(sixth 

explicit) 

Grid 

log Li 

Rate 

log L 2 

Rate 

log L 2 

Rate 

21 

-1.3472 


-1.8061 


-2.7007 


31 

-2.0866 

4.20 

-2.4296 

3.54 

-3.8229 

6.37 

41 

-2.5980 

4.09 

-2.8773 

3.58 

-4.6666 

6.75 

61 

-3.3107 

4.05 

-3.5243 

3.67 

-5.8518 

6.73 

81 

-3.8145 

4.03 

-3.9978 

3.79 

-6.6485 

6.38 


Table IH.b: Absolute error (log L 2 ) and convergence exponent with SAT parameter r = 2, for the 
fourth explicit, fourth compact and sixth explicit spatial discretizations. 

A comparison of the SAT grid refinement studies (table Ill.a and Ill.b) with those from the con- 
ventional boundary treatment (table I), indicates that the formal accuracy of the spatial operator is 
unaffected by the SAT treatment. The proof of stability given in section 1 indicated that a sufficient 
condition for stability of the scalar wave equation with the SAT method is 1 < r. The results shown 
in tables III. a and Ill.b indicate that the magnitude of the error is dependent on the value of the 
parameter r. To optimize the value of the parameter r for these simulations, the error at T = 10 
was studied as a function of r. An eigenvalue code was then used to determine the maximum CFL 
of the scheme as a function of r . The results of this study are shown in Table IV. 


r 

log L 2 

CFL 

3.0 

-3.8220 

1.17 

2.5 

-3.8221 

1.77 

2.0 

-3.8145 

2.07 

1.75 

-3.8038 

2.07 

1.50 

-3.8833 

2.07 

1.25 

-3.7460 

2.07 

1.00 

-3.6806 

2.07 

0.97 


0.0 


Table IV: Absolute error (log l 2 ) and CFL for various values of the SAT parameter r, for the fourth 
explicit spatial operator. 

Note that a fairly sharp cutoff at the theoretical value of r = 1 is observed for the fourth-order 
explicit spatial operator. (Values of r = 0.93 and r = 0.99 were obtained for the fourth-order 
compact and sixth-order explicit schemes, respectively. In addition, precise agreement was obtained 
at the r cutoff between the eigenvalue determination and the numerical simulation of the scalar wave- 
equation.) For the fourth-order explicit spatial operator, the error decreased monotonically with r. 
which suggests that the value of r should be as large as possible. Conversely, the maximum CFL 
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that is achievable with the fourth-order R-K scheme decreases dramatically at r 2. A value of 
r = 2 was determined to be optimal for these studies. 

SAT Boundary Cond itions (System) 

The last part of the validation study is to verify that the SAT boundary procedure ensures sta i ity 
for the hyperbolic system. Equation (31) defines sufficient conditions for time stability ^(2 
20 - \m\) < r < 1^(2 + 20 - |W) in terms of r and the boundary coupling matrices L 

and R. The test case chosen is the hyperbolic system 


du 

du 

1 — 

n 


dt 

dx 



dv 

dv 

0, 0 < X < 1, t > 0 

(41) 

It 

dx 



(0,0 = 

av(0, t), 

t>(l,t) = /MM)> * >° 

(42) 


u(*,0) = sin2ir*, t>(x,0) = -sin27rx, 0 < x < 1, 


The exact solution for a = /? = 1 is 


u{x, t ) 


sm2ir{x-t), v(x,t) = -sin27r(x + t), 0 < x < 1, t>0 


(44) 


The case |a 0| = 1 is neutrally stable and provides an extremely severe test of the time stability 
of a numerical method. No central difference scheme of an order greater than two, is time-stable for 
this system, in spite of the fact that the spatial operator is stable for the scalar case (a - 0 - 0). 


111 , 111 

Examples include the (3-4 3) compact and (3, 3-4-3, 3) explicit fourth-order schemes, and the ( 


6 -! 


k 2 r‘2 


ICS U1UUUC UllO -X KF J V ' ' ' 

0-0 ,0 ) sixth-order scheme that is shown in reference [4] to be time-stable for the scalar case. All 
three schemes used in the scalar analysis (fourth-order explicit and compact and sixth-order explicit), 
that satisfy the summation-by-parts property are not time-stable. In all cases, the discrete solution 
of the system defined by equations (41) through (44) diverges as time becomes large. Grid refinement 
shows Lax stability and an order property for each scheme, but not time stability. 


The scalar analysis demonstrates a precise relationship between schemes that are time-stable and 
the structure of the eigenvalue spectrum that arises from the discretization matrix. Precisely, i 
RH-P eigenvalues exist, then numerical divergence can be expected from the numerical simulation. 
Unfortunately, this statement is a function of the CFL that is used to advance the solution. (See 
reference [4].) Values of the CFL can be chosen for which no numerical divergence is experienced with 
an R-K time advancement scheme; for this reason testing the numerical stability of various spatial 

operators for the fully discrete system in time is impractical. 
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The alternative is to use the eigenvalue structure of the semi-discrete problem as the test for 
stability. If a spatial discretization operator has no RH-P eigenvalues, then it is assumed to be time- 
stable. A derivation of the discretization matrix operators for the model hyperbolic system [equations 
(41) and (42)] is presented in Appendix III. In addition, the structure of the eigenvalues is derived. 


Formir test system, we take a = 0 in (44) and thus the sufficient condition for stability becomes 

( 2— 2y/l —a 2 \ ^ ^ ^ ( 2+2\/l —a 2 \ ri' i 

1 ' - 7 - ' ^2 J* til ven a value of a and a stable scheme incorporating the SAT 

boundary treatment for the system, there exist a range in r for which the time discretization is 
stable. As in the scalar case, good agreement exists between the theoretical and numerical stability 
limit. Therefore, the agreement between the theoretical prediction and the numerical eigenvalue 
determination was used as a test of the validity of the theory. 


Table V compares the stability limits of the three high order schemes for various values of the 
parameter a; the theoretical limit is compared with that predicted from the eigenvalue determination 
for the 2x2 system. The number of grid points used in each case was 101. A study with 61 points 
showed similar results. In the study, tj is the theoretical value of r based on = T ^ an< | Tv 

is the value as determined from the eigenvalue determination. Specifically, t n was the smallest value 
of r for which the numerical eigenvalues all had negative real parts. In all cases the agreement was 
very good, which suggests the validity of the theory. 



a 

1.0 

0.99 

0.90 

0.80 

0.50 

Exact 

r T 

2.0 

1.75 

1.39 

1.25 

1.07 

fourth explicit 

tn 

2.0 

1.75 

1.39 

1.24 

1.05 

fourth compact 

tn 

2.0 

1.75 

1.39 

1.25 

1.08 

sixth explicit 

t n 

2.0 

1.72 

1.25 

1.01 

1.00 


Table V: The theoretical and numerical stability limits of SAT boundary scheme for various values 
of a. 

In these simple examples, we have demonstrated that the SAT boundary procedure retains the 
formal accuracy of the underlying spatial operator and provides a mechanism to stabilize those spatial 
operators that satisfy a summation-by-parts energy property. The resulting scheme is time-stable for 
both the scalar and system case. The numerically predicted stability boundaries for the parameter r 
closely match the theoretical predictions. From a practical perspective, the numerical stability and 
CFL of the fully discrete algorithm are functions of the value of r. The choice r = 2 seems to be 

well suited for both the scalar and system cases and guarantees stability even for the neutrally stable 
system case where a = 0 = 1 . 
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4. Conclusions 


In this paper we studied the stability and time stability of the semi-discrete hyperbolic system 
of partial differential equations. The spatial discretizations considered were high order (explicit and 
compact), and their boundary terms were constructed such that the derivative matrix satisfied a 

si unmation-by-parts formula. 

The following results were obtained: 

1. A systematic way was developed to obtain high-order accurate derivative matrices (includ- 
ing boundary terms) having a summation-by-parts property. The method is illustrated by 
finding explicit forms in the 4th order compact case. 

•2. The summation-by-parts property does not, by itself, guarantee the stability and time 
stability of the scheme, not even in the scalar case. (Refer to the explicit sixth-order 

example cited in the text.) 

3 . To overcome this difficulty we introduce the simultaneous approximation term (SAT) in 
order to account for the effect of the coupling of the physical boundary conditions. The 
SAT contains a free parameter r. 

1. We give bounds on r such that the resulting scheme for the system (or scalar) case, we 
have stability as well as time stability. 

5. Numerical studies verify the theoiy. 
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APPENDIX I 


( Const ruction of the Fourth-Order Compact Scheme 
We begin with the semi-discrete equation u t = A ] u where u = («i,u 2 , ...,u N ) , which results from a 
particular discretization of the equation u t = u x . The matrix A 1 is then decomposed as A ] = P~ X Q. 
The interior scheme used is the fourth-order compact scheme defined implicitly as 

1 du, - x + = A(u I+1 - u t -_i) (AI. 1) 

ax 4 dx 


A dx dx ' 4 dx 4 A 
Note that the interior scheme satisfies the summation-by-parts energy norm (as well as the generalized 

norm). The matrices P and Q can be written in general form, with boundary closures of arbitrary 
size iV as 



Po,o ■ 

• Po,N o 


qo,o ■ ■ 


0 ' 

p = 

Pn, o • 
0 

■ ■ PN,N \ 

I 1 I 

4 4 

; Q = 

</N, 0 • ■ 

0 

■ ■ QN,N 

-3 
4 

— 1 

col 

coirr o 


with the // matrix written as 


H = 


ho o ho. 


friV,0 

0 


N 


♦ fr/v,/v % 
x y x 


To simplify the matrix algebra, the following new matrices are introduced: 


1 0 

0 1 

-1 0 1 


C = 


4 1 0 
1 4 1 
0 1 4 





O 

o 


i 

O 

B 

L 



D = 

x y x 

0 x y x 

; A = 

0 


■ • - 


1 0 . . 0 
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Note that S, C, and D are M x M matrices, where M is an arbitrary number that corresponds to 
the number of interior points in the discretization. The structure of the matrices is tri-diagonal in 
nature. The matrix A is N x A/, and the only non-zero element is a N>1 = 1. 

Thus, we can write H , P, and Q as 


H 

x A ' 

; P = 

p 


; Q = 

Q 

w 

_xA T 

D 



C . 


-jA T 

L 4 

s 


where P, Q , and H are the N x N submatrices that involve the unknown quantities in the matrices 
P, Q , and //, respectively. 

The spatial operator that involves P and Q satisfies the generalized summation-by-parts energy 
norm if a matrix H can be found which simultaneously symmetrizes H P and yields an H Q matrix 
that is nearly skew symmetric. By defining W = H P and V = HQ, the matrices W and V become 


W = 


H P + | AA t 
x A t P + \ D A T 

4 


X AC + \H A ' 

DC + f A 7 'A _ 


V 


HQ - AA T x AS + \H .4 ' 
. x A T Q — | D A T D S + ~ A 7 A 


Thus, the matrices W and V are important to the stability properties of the spatial operator. 
Several notes about the structure of W and V should be made at this point. First, the matrices 
A A and A T A are zero except for the (N, N) and (0,0) elements, respectively. Second, the matrix 
D C + 4 A T A is automatically symmetric and it has the same tri-diagonal structure as the D and 
C matrices. Third, the matrix D & ~f A ^ A is automatically skew-symmetric which includes the 
zero at the (0,0) position. The fourth quadrant of W and V automatically satisfy the conditions on 
the generalized summation-by-parts energy norm. The remaining conditions that W and V must 
satisfy, written in terms of the submatrices H, P, Q, C, D, .S’, and A, are 

HP = (H P) T 


(Al. 2) 



(AI. 3) 


-A t H t + xC t A t = \d A t + x A t P 
4 4 

H Q + {H Qf = ^ A A t + A 8 0 , 0 I (AI. 4) 

-A t H t + xS t A t = '-DA t - xA t Q (AI. 5) 

4 4 

where A 8 0 ,o is the non-zero element that occurs in the first row and column of the matrix. This 
contribution to equation (AI. 4 ) allows for a non-zero value at the (0,0) element in the matrix V . 

By expanding the specific terms in equations (AI. 2 through AI. 5), we have 



ho, n 

. • hn,n 


Pn, 0 

Pn,n 


0 

0 


0 

0 

A T H T = 

0 

0 

II 

L 0 

0 


A T Q = 

9n,0 

0 

qn,n 

0 

II 

hi 

h 

' 0 . 

. 01 ' 

1 

4 

0 


0 

0 


0 

0 _ 


' 0 . 

. 00 ' 


' 0 . 

1 

s* 

o 


3 



X 


4 

0 

; DA T = 


0 

0 

0 _ 


0 

0 


By comparing the matrices involved in equation (AI. 3), it is apparent that 

7 h k ,N + x8 k ,N = x p N ,k + 7 h,N, k = 0,N (AI. 6 ) 

4 4 

Similarly, equation (AI. 5) yields the expression 

7 h kl N = — x qN k + ~t~ &k,N] k = 0, N (AI. 7) 

4 4 

Eliminating between equation (AI. 6 ) and equation (AI. 7) yields the expression 

9tv, it = —3 pTv.fc + 3 8 k , n] k = 0,N 

24 


(AI- 8 ) 



These properties of the matrices P and Q must be satisfied regardless of the order properties of the 
boundary. 

We now derive the additional constraints that must he satisfied near the boundaries to guarantee 
the order properties of these points. Substitution of the equations Uj = j T and ^ = r j r ~ x into 
the matrices Q and f\ respectively, yields the constraints that ensure the accuracy of the boundary 
points. The general expression at the boundary written in terms of an arbitrary accuracy r becomes 


;V 

r E Pk,j J 

7 = 0 


r-t 


+ - f>h,N (N + 1 ) 


r — 1 


E 'A-., / + U,, n (N+\) t : k = 0..V 

7=0 ' 1 


(Al. 9) 


Third-order accuracy at the boundary points requires r = 0,3 with N > 3. 


Thus far, we have not specified the exact value of the parameter N . We now specify a precise 
value for the parameter N so that specific boundary conditions can be derived for the fourth-order 
interior Pade scheme. To retain the formal accuracy of the interior scheme, the boundary closure 
must be accomplished to at least third-order accuracy, and requires that N > 3. For N — 3. 
equation (AI. 9) can be written concisely in matrix notation as 


P 


' 0 + 0" 1 

i 

* 0° 

2 * 0 1 

3 * 0 2 ' 


0° 

O 1 

0 2 

0 ;} 1 


■ 0 

0 

0 

0 ■ 

0 * l- 1 

i 

* i° 

2 * l 1 

3 * l 2 


i° 

l 1 

l 2 

1 } 


0 

0 

0 

0 

0 * 2 -1 

i 

* 2° 

2 * 2 1 

3 * 2 2 

= Q 

2° 

2 1 

2 2 

2 3 

+ 

0 

0 

0 

0 

. o * 

i 

*3° 

2 * 3 1 

3 * 3 2 . 


3° 

3 1 

3 2 

3 3 


3 

- 1 

11 

4 

10 

3(> . 


Solving this expression for the matrix Q results in the expression 


Q = P 



+ 


0 0 0 0 ■ 

0 0 0 0 

0 0 0 0 

7 -5 17 -23 

24 4 8 12 - 


which relates the matrix Q to the matrix P through third-order accuracy constraints. 

We will now solve for the last row of the matrices P and Q and for the last column of the matrix 
//. Equation (Ah 9) is written for k — N, and and (defined in equation (AI. S)) arc used 
to yield the relationship 



r = 0,3 (AI. 10 ) 


r 51 PNj f 1 + 7 ( N + !) r 1 = -3 E j r + T (W + l) r + 3 N r ; 

j=o 4 j=o 


Setting /V = 3 and solving the system for p 3) fc, k = 0,3 yields p 3?0 = P3,i = 0 , P3,2 = 4? an ^ p 3>3 = 1 . 
Equation (AI. 8 ) can be used to show that q 3i0 = 93,1 = 0 , q 3 ;i = — §> an d 93,3 = 0 . Similarly, 
equation (AI. 6) yields the values of h *^3 as /io,3 = ^1,3 — 0 , ^2,3 = an d /i 3)3 = y. Thus, the last 
row of P and Q are the same as the interior scheme. In addition, the specific form of the matrix H 
must be 


^0,2 0 

^1,2 0 

^2,2 ^ 

^3,2 y . 

Thus, accuracy constraints on the last row of the matrices P and Q, combined with the structure 
requirements imposed by equations (AI. 3) and (AI. 5 ), allow for the direct solution of the last rows 
of P and Q , and the last column of H . Multiplying the expression relating P to Q by the matrix H , 
and using the substitutions H Q = V and H P = W yields the expression for V of the form 


H 


ho } o ^0,1 
^1,0 ^1,1 
^2,0 ^2,1 
^ 3,0 ^ 3,1 



r -11 
6 

3 

-3 

2 

1 -1 

3 


■ 0 

0 

0 

0 ‘ 


-1 

-1 

1 

-1 


0 

0 

0 

0 

li 

3 

1 

2 

1 

6 

1 

+ 

7 x 

—5 x 

17 x 

-23 x 


6 

-1 

L 3 

— I 

3 

2 

2 

—3 

1 

'Iconics 


l^ls 

1 

4 

-$y 

4 

i±L 

8 

-tty 
12 J 


Solving for W and V such that equation (AI. 2) (where W — W T ) and equation (AI. 4 ) (V + V — 


A T + . 

A £o,o /) are satisfied 

to obtain 




—9 a 

1536^+1536/3—899 a 

768 7+768/3-703 a 

15367+1536/3-1481 a 


16 

768 

192 

768 


1 536 -7+1536 0—899 a 

0 

15367+1536 /3-1277a 

768 7+768/3-733 a 


768 

256 

192 

V 7 = 

768 7 + 768/3-703 a 

1 536 7 + 1536/3- 1277 a 

0 

15367 + 1536 /3-947a 


192 

256 

768 


1536 7+1536/3-1481 a 

768 -y+768 /?— 733 a 

15367+1536/3-947 a 

-3 a 


L 768 

192 

768 

32 


and 
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137 a — 192 7 

192 


— 512 '> — 128/3+525 or 
1 28 


145 o — 144 7 
48 


w 


-5127-128/3+5250 

128 

145 0—144 7 

48 


-816 7-384/3 + 737 0 
48 


-3072 7-1344 /3+3Q01 o 
192 


-3072 7-1344 /3+3001 o -3264 7-153 6 /j+3011 a 


192 


192 


557 q — 576 7 

T02 

-15367-384 ,/3+1543o 
384 


0 


557o-576 7 
192 


-15367-384 /3+ 1543 o 
384 


with x — and y = a. Three arbitrary parameters remain after all accuracy, symmetry and 
skew-symmetry conditions are satisfied. 

The final step in the discretization is to find a specific form of the matrix P that will lead to a 
simple algorithm. Because the matrix P is tri-diagonal in the interior, the boundary closure should 
retain the tri-diagonal structure. After P is specified, we can solve for the matrix H from // = V / J_1 
if the inverse of P exists, and the last column of H is [0,0,y,^] r . The matrix Q follows immediately 
from Q — P V~ } W . The last test is to ensure that both W and that the full matrix 14 are positive' 
definite. 

Many matrices P have been found that satisfy all of the criteria given in the generalized summation- 
by-parts energy norm analysis. From a numerical perspective, all behaved similarly. The results 
presented here are those that were the simplest to code. Choosing a specific matrix P of the form 


r in 

429 


0 0 


P 


3563 -1 n 

1688 8 U 

43 1893 139 

17 1054 186 


0 0 i 1 

4 


yields a matrix Q of the form 


-289 

279 

75 

— 7 

234 

286 

286 

2574 

-8635 

6987 

1851 

-203 

3376 

3376 

3376 

3376 

-15043 

-4089 

147 

29353 

18972 

2108 

124 

18972 

0 

0 

-3 

A 

0 


The resulting matrix H is therefore 
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70282007653 

-9426299 

-192913 

0 

7658388480 

2268480 

1067520 

-55530689643 

8051589 

149823 

0 

2552796160 

756160 

355840 

63842626133 

-9153739 

-4433 

-1 

2552796160 

7561 60 

355840 

8 

-71498870443 

10110149 

102703 

1 

7658388480 

2268480 

1067520 

J 


^Frorn a practical point of view, the inconvenient form of the // matrix is not of great concern 
since the matrix H is only inverted once and one column is stored for use. 

The matrices P, Q, and H can he used to establish both the symmetry of the matrix V and the 
near skew-symmetry of the matrix W . The first six rows and columns of the V matrix are 


16513 

-261 

2993 

-6223 

0 

0 

46080 

5120 

15360 

46080 

-261 

9153 

-2943 

1611 

0 

0 

5120 

5120 

5120 

5120 

2993 

-2943 

7473 

-2063 

-1 

0 

15360 

5120 

5120 

1 5360 

32 

-6223 

1611 

-2063 

47953 

1 

-1 

4 ( i 080 

5120 

1 5360 

46080 

8 

32 

0 

0 

-1 

1 

15 

1 

32 

8 

16 

8 

0 

0 

0 

-1 

1 

75 

32 

8 

16 


The first six rows and columns of the W matrix are 



- ^9 

45 

-11 

-7 

0 

0 


16 

64 

128 

128 


-45 

64 

0 

81 

128 

9 

128 

0 

0 


11 

-81 

0 

41 

-3 

0 


128 

128 

64 

32 

W = 

7 

-9 

-41 

0 

3 

-3 


128 

128 

64 

4 

32 


0 

0 

3 

32 

-3 

4 

0 

3 

4 


0 

0 

0 

3 

32 

^3 

4 

0 


As shown, the matrix W is nearly skew symmetric, and the matrix V is symmetric. For the matrix W 
is positive definite, it is necessary to show that every submatrix is positive definite. The inner scheme 
is diagonally dominant and contributes to the definiteness of the complete matrix W . However, the 
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boundary elements are not diagonally dominant, and suppress the positive-definiteness. The 1 X I 
boundary matrix W' = It + | A A 1 has the following characteristic polynomial 

751971720 A 1 - 2507814400 A 1 + 5299008028 A’ 2 - 2079222121 A + 520779791 = 0 (Al. 11) 

The symmetry of the W f matrix and the alternating signs of each term in the characteristic polynomial 
guarantee that the matrix is positive definite. 1 lie characteristic polynomial of every submatrix (up 
to ten points, which includes four boundary and six interior points) of the matrix It results in a 
positive definite matrix. No proof that the complete discretization is positive definite for an arbitrary 
number of interior points lias been found. 

The accuracy of the new scheme is third order at the boundaries and lourth order in the interior. 
To show this, the Taylor expansion tor long wavelength modes is made using the stencil at each of 
the first, four points. The results are 


*4 + 


17 

640 


it + 

»£ - 


42 


2016 

98 


2005 


y - 


1 180 


■c + • • • 
y + . . . 
y + . . . 
y +... 


(AI. 12) 


At high resolution, the boundary points behave with third-order truncation error; the interior behaves 
with fourth-order error. Therefore, the resulting scheme is formally fourth-order accurate'. 



APPENDIX II 


Sixth Order Explicit Scheme 

Here, we derive an explicit scheme that is formally sixth-order accurate. Unlike the fourth-order 
compact case presented earlier, the matrix H can be the identity matrix. To constrain the matrix P 
to be symmetric and the matrix Q be nearly skew symmetric, six alternative formulas are required 
at the boundaries, each of which is closed to fifth-order accuracy to retain the formal accuracy. The 
corner 7x7 submatrices of the global matrices P and Q can be written as 


Qe = 



2113 

18487 

553 

14759 

(-29269) 

54839 

0 " 



10800 

345600 

57600 

172800 

1 72800 

345600 



18487 

345600 

175781 

51840 

(-28361) 

6912 

129329 

34560 

(-346319) 

207360 

(-19061) 

172800 

0 



553 

(-28361) 

43807 

(-915) 

126833 

(-39307) 

0 



57600 

6912 

5184 

128 

34560 

518400 



14759 

172800 

1 29329 
34560 

(~ 915 ) 

128 

67769 

8640 

(-25289) 

6912 

34811 

172800 

0 



(-29269) 
1 72800 

(-346319) 

207360 

1 26833 
34560 

(-25289) 

6912 

156053 

51840 

(-21059) 

115200 

0 



54839 

345600 

(-19061) 

172800 

(-39307) 

518400 

34811 

172800 

(-21059) 

115200 

32569 

32400 

0 



0 

0 

0 

0 

0 

0 

1 



Lzll 

2 

1235503 

1036800 

(-859597) 

518400 

398 

225 

(-603059) 

518400 

14969 

41472 

0 

t 

-1235503) 

1036800 

0 

16343 

5760 

(-68005) 

20736 

186797 

69120 

(-184657) 

172800 

0 


859597 

518400 

(-16343) 

5760 

0 

128759 

51840 

(-18743) 

6912 

3799 

2700 

0 


(-398) 

225 

68005 

20736 

(-128759) 

51840 

0 

110351 

51840 

(-607693) 

518400 

1 

60 


603059 

518400 

(-186797) 

69120 

18743 

6912 

(-110351) 

51840 

0 

376549 

345600 

iz2l 

20 

( 

-14969) 

41472 

184657 
1 72800 

(-3799) 

2700 

607693 

518400 

(-376549) 

345600 

0 


3 

4 


0 

0 

0 

lzll 

60 

3 

20 

iz 31 

4 

0 


die characteristic polynomial of the matrix P& is 

10399739562845798400000000 A 6 - 248512609916244983808000000 A 5 
+ 1003578630643249838161920000 A 4 - 1639038223377237368051712000 A 3 
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+ 1248376737213799711434406800 A' 2 - 412235365042816633559197440 A 

+ 37455444120716264727507839 = 0 (AIL 1) 

The symmetry of the matrix P 6 and the alternating signs of the terms in the polynomial are sufficient 
for positive definiteness of both the matrix P 6 and the global matrix P. 


The truncation error at the boundary points is 


6448299997451547397244 467^ 

^ + 224732664724297588365047034 

55 17 8459341 997062554732 1£ 6 
^ " 1 1 2366332362148794 1 825235 1 70 

903781 1404281 6098962729 6 19£ 6 
~ 2247326647242975883650470340 

62520732887440126777806839^ 
^ + 2247326647242975883650470340 

21521021082694965917733 H 6 
^ + 1123663323621487941825235170 

71015802541971 16302053905^ 

^ ~ 224732664724297588365047034 


(All. 2) 


which indicates fifth-order accuracy at the six boundary points. 
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APPENDIX III 


Eigenvalues of the Discrete System 

Hie eigenvalues of the semi-discrete system are used in the results section to compare the theoretical 
and the numerical stability boundaries. The model equation is the hyperbolic system used in the 
main text and defined by equations (41), and (42). For convenience, we define the' (N + 1) x(N + l) 

mat lix A = P 1 Q. I he matrix A contains all the information from the spatial discretization 
operator . I lu semi-discrete form of equation (41) becomes 

du 

It + Au = °’ 

dv 

di ~ Av = °> (Ain. i) 

with the boundary conditions defined by equation (42). In matrix notation, the discrete system takes 
the form 


du 

Ot 


A* . aB 
8J~ l BJ . J-'sVj 


where 


til 


a 1,1 

a l,2 

■ • a l,N - 1 

a \,N 

UN- 1 


a 2,1 

«2,2 

■ a 2,N-l 

U'2 y N 

Ujsj 

Vq 

; A* = 





V\ 


«JV-1,1 

a N-\;2 

■ ■ «AT-l,A r -l 

a N-\,N 

. ^V-l 


a N,\ 

a N, 2 * - 

' • d-N,N - 1 

a N,N 


and 


B = 


°I,0 

0 

0 . 

..o' 


- 

1 ' 

«2,0 

0 

0 . 

. 0 


0 

1 





; J = 



«jv- i,o 

0 

0 . . 

. 0 


I 

0 

a NM 

0 

0 . . 

. 0 


I 
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Note that JJ = 1, so that J = J~\ The vector u is the concatenated vector of discrete values 
from the scalar vectors u and v with the elements u 0 and i>n removed. These elements are removed 
because the physical boundary condition relates them to known elements in the vector u, so that and 
need not need to solve for them. The matrix A f is the N x N submatrix of A which is obtained by 
eliminating the zeroth row and zeroth column. Note that this was the matrix that was analyzed in the 
scalar analysis to determine time stability of the spatial operator. The matrix B is zero everywhere 
except the first column, where the zeroth column of the original A matrix is written. This column is 
precisely the coupling between the u and v vector, which occurs at the boundary. 

It is instructive to relate the system eigenvalues to those obtained in the scalar analysis [(/l 1 
XI) u = 0]. By defining the matrix H~' and H as 

" y/fil . y/aJ 1 [ • -V" 1 

H -» = -J— • • • ; H = 7= • ' a 

[ -ytfl ■ sfeJ J V 2 V^ L VPJ • '/$ J 

with // -1 H — H = /, we note that the system matrix can be made block diagonal with the 

similarity transform H 

y/fil . s/aJ "1 T /1 + . otB y/al . -yftil 

-VPI ■ J L PJ-'BJ ■ J-M+J J [y/PJ ■ y/PJ 

A' + y/^BJ . 0 

0 . A* - yf^BJ 

For scalar time-stable spatial schemes, the eigenvalues of the matrix A j are bounded to the left half- 
plane. Note that for a = 0 (or f) = 0) the contribution from the boundary coupling matrix B is 
identically zero, and the eigenvalues of the resulting system are simply the scalar eigenvalues with 
a multiplicity of two. For non-zero values of the parameters cv and /?, the eigenvalues of the total 
matrix are different from those of the original matrix A f . Also note that two distinct eigenvalue 
scenarios exist for the boundary parameters a and /?, depending on whether their signs are equal or 
opposite. 
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